// super_math.h
#ifndef SUPER_MATH_H
#define SUPER_MATH_H

#include <cmath>
#include <complex>
#include <vector>
#include <functional>
#include <algorithm>
#include <numeric>
#include <random>
#include <type_traits>
#include <concepts>
#include <numbers>
#include <iostream>
#include <iomanip>
#include <limits>

using namespace std;

// ==================== 类型定义和概念 ====================

template<typename T>
concept Arithmetic = is_arithmetic_v<T>;

template<typename T>
concept FloatingPoint = is_floating_point_v<T>;

template<typename T>
concept Complex = is_same_v<T, complex<float>> || 
                  is_same_v<T, complex<double>> || 
                  is_same_v<T, complex<long double>>;

// 高精度实数类型
using BigFloat = long double;

// 数学常量
namespace MathConstants {
    constexpr double PI = 3.14159265358979323846;
    constexpr double E = 2.71828182845904523536;
    constexpr double GOLDEN_RATIO = 1.61803398874989484820;
    constexpr double EULER_MASCHERONI = 0.57721566490153286060;
    constexpr double CATALAN = 0.91596559417721901505;
    constexpr double GLAISHER_KINKELIN = 1.28242712910062263687;
    
    // 物理常数
    constexpr double SPEED_OF_LIGHT = 299792458.0;
    constexpr double GRAVITATIONAL_CONSTANT = 6.67430e-11;
    constexpr double PLANCK_CONSTANT = 6.62607015e-34;
    constexpr double BOLTZMANN_CONSTANT = 1.380649e-23;
    constexpr double AVOGADRO_CONSTANT = 6.02214076e23;
}

// ==================== 基础数学函数扩展 ====================

namespace SuperMath {
    
    // ==================== 特殊函数 ====================
    
    // Gamma函数
    template<FloatingPoint T>
    T gamma(T x) {
        // Lanczos近似
        constexpr T g = 7.0;
        constexpr T coefficients[] = {
            0.99999999999980993, 676.5203681218851, -1259.1392167224028,
            771.32342877765313, -176.61502916214059, 12.507343278686905,
            -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
        };
        
        if (x < 0.5) {
            return PI / (sin(PI * x) * gamma(1.0 - x));
        }
        
        x -= 1.0;
        T result = coefficients[0];
        for (int i = 1; i < sizeof(coefficients)/sizeof(coefficients[0]); ++i) {
            result += coefficients[i] / (x + T(i));
        }
        T t = x + g + 0.5;
        return sqrt(2.0 * PI) * pow(t, x + 0.5) * exp(-t) * result;
    }
    
    // Beta函数
    template<FloatingPoint T>
    T beta(T x, T y) {
        return gamma(x) * gamma(y) / gamma(x + y);
    }
    
    // 误差函数
    template<FloatingPoint T>
    T erf(T x) {
        // 近似计算
        T t = 1.0 / (1.0 + 0.3275911 * abs(x));
        T poly = t * (0.254829592 + t * (-0.284496736 + t * (1.421413741 + 
                t * (-1.453152027 + t * 1.061405429))));
        T result = 1.0 - poly * exp(-x * x);
        return (x >= 0) ? result : -result;
    }
    
    // 互补误差函数
    template<FloatingPoint T>
    T erfc(T x) {
        return 1.0 - erf(x);
    }
    
    // Bessel函数 (第一类)
    template<FloatingPoint T>
    T bessel_j(int n, T x) {
        if (x == 0.0) return (n == 0) ? 1.0 : 0.0;
        
        const int max_iter = 100;
        T sum = 0.0;
        T term = 1.0;
        T x_half = x / 2.0;
        
        for (int k = 0; k < max_iter; ++k) {
            if (k >= n) {
                sum += term;
            }
            term *= -x_half * x_half / ((k + 1) * (k + n + 1));
            if (abs(term) < 1e-15) break;
        }
        
        return pow(x_half, n) * sum / tgamma(n + 1.0);
    }
    
    // ==================== 数论函数 ====================
    
    // 判断素数
    template<Arithmetic T>
    bool is_prime(T n) {
        if (n < 2) return false;
        if (n == 2 || n == 3) return true;
        if (n % 2 == 0 || n % 3 == 0) return false;
        
        for (T i = 5; i * i <= n; i += 6) {
            if (n % i == 0 || n % (i + 2) == 0) return false;
        }
        return true;
    }
    
    // 欧拉函数
    template<Arithmetic T>
    T euler_phi(T n) {
        T result = n;
        for (T i = 2; i * i <= n; ++i) {
            if (n % i == 0) {
                while (n % i == 0) n /= i;
                result -= result / i;
            }
        }
        if (n > 1) result -= result / n;
        return result;
    }
    
    // 莫比乌斯函数
    template<Arithmetic T>
    int mobius(T n) {
        if (n == 1) return 1;
        
        T temp = n;
        int prime_count = 0;
        
        for (T i = 2; i * i <= temp; ++i) {
            if (temp % i == 0) {
                ++prime_count;
                temp /= i;
                if (temp % i == 0) return 0;
            }
        }
        if (temp > 1) ++prime_count;
        
        return (prime_count % 2 == 0) ? 1 : -1;
    }
    
    // 因数个数
    template<Arithmetic T>
    T divisor_count(T n) {
        T count = 1;
        for (T i = 2; i * i <= n; ++i) {
            T power = 0;
            while (n % i == 0) {
                n /= i;
                ++power;
            }
            if (power > 0) count *= (power + 1);
        }
        if (n > 1) count *= 2;
        return count;
    }
    
    // 因数之和
    template<Arithmetic T>
    T divisor_sum(T n) {
        T sum = 1;
        for (T i = 2; i * i <= n; ++i) {
            T power = 0;
            T term = 1;
            while (n % i == 0) {
                n /= i;
                ++power;
                term = term * i + 1;
            }
            if (power > 0) sum *= term;
        }
        if (n > 1) sum *= (1 + n);
        return sum;
    }
    
    // ==================== 组合数学 ====================
    
    // 阶乘
    template<Arithmetic T>
    T factorial(T n) {
        if (n < 0) return 0;
        T result = 1;
        for (T i = 2; i <= n; ++i) {
            result *= i;
        }
        return result;
    }
    
    // 双阶乘
    template<Arithmetic T>
    T double_factorial(T n) {
        if (n <= 1) return 1;
        T result = 1;
        for (T i = n; i > 1; i -= 2) {
            result *= i;
        }
        return result;
    }
    
    // 组合数
    template<Arithmetic T>
    T nCr(T n, T r) {
        if (r < 0 || r > n) return 0;
        if (r > n - r) r = n - r;
        
        T result = 1;
        for (T i = 1; i <= r; ++i) {
            result = result * (n - r + i) / i;
        }
        return result;
    }
    
    // 排列数
    template<Arithmetic T>
    T nPr(T n, T r) {
        if (r < 0 || r > n) return 0;
        T result = 1;
        for (T i = 0; i < r; ++i) {
            result *= (n - i);
        }
        return result;
    }
    
    // 卡特兰数
    template<Arithmetic T>
    T catalan(T n) {
        return nCr(2 * n, n) / (n + 1);
    }
    
    // ==================== 数值分析 ====================
    
    // 数值积分 - 辛普森法则
    template<FloatingPoint T, typename Func>
    T integrate_simpson(Func f, T a, T b, int n = 1000) {
        T h = (b - a) / n;
        T sum = f(a) + f(b);
        
        for (int i = 1; i < n; i += 2) {
            sum += 4.0 * f(a + i * h);
        }
        for (int i = 2; i < n; i += 2) {
            sum += 2.0 * f(a + i * h);
        }
        
        return sum * h / 3.0;
    }
    
    // 数值积分 - 高斯求积
    template<FloatingPoint T, typename Func>
    T integrate_gauss(Func f, T a, T b, int n = 10) {
        // 高斯点权重 (n=10)
        constexpr T points[] = {
            -0.9739065285171717, -0.8650633666889845, -0.6794095682990244,
            -0.4333953941292472, -0.1488743389816312, 0.1488743389816312,
            0.4333953941292472, 0.6794095682990244, 0.8650633666889845,
            0.9739065285171717
        };
        constexpr T weights[] = {
            0.0666713443086881, 0.1494513491505806, 0.2190863625159820,
            0.2692667193099963, 0.2955242247147529, 0.2955242247147529,
            0.2692667193099963, 0.2190863625159820, 0.1494513491505806,
            0.0666713443086881
        };
        
        T sum = 0.0;
        T scale = (b - a) / 2.0;
        T shift = (a + b) / 2.0;
        
        for (int i = 0; i < n; ++i) {
            sum += weights[i] * f(scale * points[i] + shift);
        }
        
        return sum * scale;
    }
    
    // 数值微分
    template<FloatingPoint T, typename Func>
    T derivative(Func f, T x, T h = 1e-8) {
        return (f(x + h) - f(x - h)) / (2.0 * h);
    }
    
    // 二阶数值微分
    template<FloatingPoint T, typename Func>
    T second_derivative(Func f, T x, T h = 1e-6) {
        return (f(x + h) - 2.0 * f(x) + f(x - h)) / (h * h);
    }
    
    // 牛顿法求根
    template<FloatingPoint T, typename Func, typename Deriv>
    T newton_method(Func f, Deriv df, T x0, int max_iter = 100, T tol = 1e-12) {
        T x = x0;
        for (int i = 0; i < max_iter; ++i) {
            T fx = f(x);
            T dfx = df(x);
            
            if (abs(dfx) < tol) break;
            
            T delta = fx / dfx;
            x -= delta;
            
            if (abs(delta) < tol) break;
        }
        return x;
    }
    
    // ==================== 线性代数 ====================
    
    // 矩阵乘法
    template<Arithmetic T>
    vector<vector<T>> matrix_multiply(const vector<vector<T>>& A, 
                                     const vector<vector<T>>& B) {
        int n = A.size(), m = B[0].size(), p = B.size();
        vector<vector<T>> C(n, vector<T>(m, 0));
        
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                for (int k = 0; k < p; ++k) {
                    C[i][j] += A[i][k] * B[k][j];
                }
            }
        }
        return C;
    }
    
    // 矩阵行列式
    template<Arithmetic T>
    T determinant(vector<vector<T>> matrix) {
        int n = matrix.size();
        T det = 1;
        
        for (int i = 0; i < n; ++i) {
            // 寻找主元
            int pivot = i;
            for (int j = i + 1; j < n; ++j) {
                if (abs(matrix[j][i]) > abs(matrix[pivot][i])) {
                    pivot = j;
                }
            }
            
            if (pivot != i) {
                swap(matrix[i], matrix[pivot]);
                det = -det;
            }
            
            if (abs(matrix[i][i]) < 1e-12) return 0;
            
            det *= matrix[i][i];
            
            // 消元
            for (int j = i + 1; j < n; ++j) {
                T factor = matrix[j][i] / matrix[i][i];
                for (int k = i + 1; k < n; ++k) {
                    matrix[j][k] -= factor * matrix[i][k];
                }
            }
        }
        
        return det;
    }
    
    // 矩阵转置
    template<Arithmetic T>
    vector<vector<T>> transpose(const vector<vector<T>>& matrix) {
        int n = matrix.size(), m = matrix[0].size();
        vector<vector<T>> result(m, vector<T>(n));
        
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < m; ++j) {
                result[j][i] = matrix[i][j];
            }
        }
        return result;
    }
    
    // ==================== 统计函数 ====================
    
    // 平均值
    template<Arithmetic T>
    T mean(const vector<T>& data) {
        return accumulate(data.begin(), data.end(), T(0)) / data.size();
    }
    
    // 方差
    template<FloatingPoint T>
    T variance(const vector<T>& data, bool population = true) {
        T m = mean(data);
        T sum_sq = 0.0;
        for (const auto& x : data) {
            sum_sq += (x - m) * (x - m);
        }
        return sum_sq / (data.size() - (population ? 0 : 1));
    }
    
    // 标准差
    template<FloatingPoint T>
    T standard_deviation(const vector<T>& data, bool population = true) {
        return sqrt(variance(data, population));
    }
    
    // 相关系数
    template<FloatingPoint T>
    T correlation(const vector<T>& x, const vector<T>& y) {
        T mean_x = mean(x), mean_y = mean(y);
        T numerator = 0.0, denom_x = 0.0, denom_y = 0.0;
        
        for (size_t i = 0; i < x.size(); ++i) {
            numerator += (x[i] - mean_x) * (y[i] - mean_y);
            denom_x += (x[i] - mean_x) * (x[i] - mean_x);
            denom_y += (y[i] - mean_y) * (y[i] - mean_y);
        }
        
        return numerator / sqrt(denom_x * denom_y);
    }
    
    // 线性回归
    template<FloatingPoint T>
    pair<T, T> linear_regression(const vector<T>& x, const vector<T>& y) {
        T mean_x = mean(x), mean_y = mean(y);
        T numerator = 0.0, denominator = 0.0;
        
        for (size_t i = 0; i < x.size(); ++i) {
            numerator += (x[i] - mean_x) * (y[i] - mean_y);
            denominator += (x[i] - mean_x) * (x[i] - mean_x);
        }
        
        T slope = numerator / denominator;
        T intercept = mean_y - slope * mean_x;
        
        return {slope, intercept};
    }
    
    // ==================== 信号处理 ====================
    
    // 离散傅里叶变换
    template<FloatingPoint T>
    vector<complex<T>> dft(const vector<T>& signal) {
        int N = signal.size();
        vector<complex<T>> result(N);
        
        for (int k = 0; k < N; ++k) {
            complex<T> sum(0, 0);
            for (int n = 0; n < N; ++n) {
                T angle = -2.0 * PI * k * n / N;
                sum += signal[n] * complex<T>(cos(angle), sin(angle));
            }
            result[k] = sum;
        }
        return result;
    }
    
    // 快速傅里叶变换
    template<FloatingPoint T>
    vector<complex<T>> fft(const vector<complex<T>>& signal) {
        int N = signal.size();
        if (N <= 1) return signal;
        
        // 分割奇偶索引
        vector<complex<T>> even(N/2), odd(N/2);
        for (int i = 0; i < N/2; ++i) {
            even[i] = signal[2*i];
            odd[i] = signal[2*i + 1];
        }
        
        // 递归计算
        vector<complex<T>> even_fft = fft(even);
        vector<complex<T>> odd_fft = fft(odd);
        
        // 合并结果
        vector<complex<T>> result(N);
        for (int k = 0; k < N/2; ++k) {
            T angle = -2.0 * PI * k / N;
            complex<T> twiddle = complex<T>(cos(angle), sin(angle)) * odd_fft[k];
            result[k] = even_fft[k] + twiddle;
            result[k + N/2] = even_fft[k] - twiddle;
        }
        
        return result;
    }
    
    // ==================== 插值函数 ====================
    
    // 线性插值
    template<FloatingPoint T>
    T linear_interpolate(T x, T x0, T y0, T x1, T y1) {
        return y0 + (y1 - y0) * (x - x0) / (x1 - x0);
    }
    
    // 拉格朗日插值
    template<FloatingPoint T>
    T lagrange_interpolate(T x, const vector<T>& x_points, const vector<T>& y_points) {
        T result = 0.0;
        int n = x_points.size();
        
        for (int i = 0; i < n; ++i) {
            T term = y_points[i];
            for (int j = 0; j < n; ++j) {
                if (j != i) {
                    term *= (x - x_points[j]) / (x_points[i] - x_points[j]);
                }
            }
            result += term;
        }
        return result;
    }
    
    // ==================== 几何函数 ====================
    
    // 点到直线距离
    template<FloatingPoint T>
    T point_to_line_distance(T x, T y, T x1, T y1, T x2, T y2) {
        T A = y2 - y1;
        T B = x1 - x2;
        T C = x2 * y1 - x1 * y2;
        return abs(A * x + B * y + C) / sqrt(A * A + B * B);
    }
    
    // 多边形面积
    template<FloatingPoint T>
    T polygon_area(const vector<pair<T, T>>& points) {
        T area = 0.0;
        int n = points.size();
        
        for (int i = 0; i < n; ++i) {
            int j = (i + 1) % n;
            area += points[i].first * points[j].second;
            area -= points[j].first * points[i].second;
        }
        
        return abs(area) / 2.0;
    }
    
    // ==================== 随机分布 ====================
    
    // 正态分布
    template<FloatingPoint T>
    T normal_pdf(T x, T mu = 0.0, T sigma = 1.0) {
        return exp(-0.5 * pow((x - mu) / sigma, 2)) / (sigma * sqrt(2.0 * PI));
    }
    
    // 正态分布CDF
    template<FloatingPoint T>
    T normal_cdf(T x, T mu = 0.0, T sigma = 1.0) {
        return 0.5 * (1.0 + erf((x - mu) / (sigma * sqrt(2.0))));
    }
    
    // 泊松分布
    template<Arithmetic T>
    T poisson_pmf(T k, T lambda) {
        return exp(-lambda) * pow(lambda, k) / factorial(k);
    }
    
    // 指数分布
    template<FloatingPoint T>
    T exponential_pdf(T x, T lambda) {
        return lambda * exp(-lambda * x);
    }
    
} // namespace SuperMath

// ==================== 便捷函数宏 ====================

#define SMATH SuperMath

// 数学常量快捷访问
constexpr double PI = MathConstants::PI;
constexpr double E = MathConstants::E;
constexpr double GOLDEN_RATIO = MathConstants::GOLDEN_RATIO;

#endif // SUPER_MATH_H